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Abstract 


Using archival Fermi-LAT data with a time span of ~12 yr, we study the population of Millisecond Pulsars 
(MSPs) in Globular Clusters (GICs) and investigate their dependence on cluster dynamical evolution in the Milky 
Way. We show that the y-ray luminosity (L,) and emissivity (i.e., є; = L,/M, with M the cluster mass) are good 
indicators of the population and abundance of MSPs in GICs, and they are highly dependent on the dynamical 
evolution history of the host clusters. Specifically speaking, the dynamically older GICs with more compact 
structures are more likely to have larger L., and €,, and these trends can be summarized as strong correlations with 
cluster stellar encounter rate Г and the specific encounter rate (A = Г/М), with L, ox [409 and €, ex A0730 
for dynamically normal GICs. However, as GICs evolve into deep core collapse, these trends are found to be 
reversed, implying that strong encounters may have lead to the disruption of Low-Mass X-ray Binaries and 
ejection of MSPs from core-collapsed systems. Besides, the GICs are found to exhibit larger e, with increasing 
stellar mass function slope (є, cx ПО es decreasing tidal radius (e, cx M and distances from the 


Galactic Center (GC, є « Ro” +0.21). These correlations indicate that, as GICs losing kinetic energy and spiral in 
toward the GC, tidal stripping and mass segregation have a preference in leading to the loss of normal stars from 
GICs, while MSPs are more likely to concentrate to cluster center and be deposited into the GC. Moreover, we 
gauge e, of GICs is ~10-1000 times larger than the Galactic bulge, the latter is thought to reside thousands of 
unresolved MSPs and may be responsible for the GC ^-ray excess, which supports that GICs are generous 
contributors to the population of MSPs in the GC. 
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A Fermi-LAT Study of Globular Cluster Dynamical Evolution in the Milky 


1. Introduction 


Globular Clusters (GICs) are self-gravitating systems that evolve 
with dynamical timescale (1.е., two-body relaxation timescale) 
much smaller than their age (Heggie & Hut 2003), and the binary 
burning process (ie. extraction of binary binding energy via 
dynamical encounters) is the “heating mechanism" that supports 
the cluster against gravother-thermal collapse (Fregeau et al. 2003), 
which also makes GICs promising breeding grounds for exotic 
objects, such as blue straggler stars (BSS, Fregeau et al. 2004; 
Chatterjee et al. 2013), coronal active binaries (ABs), cataclysmic 
variables (CVs; Ivanova et al. 2006; Shara & Hurley 2006; Belloni 
et al. 2016, 2017, 2019), low-mass X-ray binaries (LMXBs; Rasio 
et al. 2000; Ivanova et al. 2008; Kremer et al. 2018), millisecond 
pulsars (MSPs; the offspring of LMXBs; Ye et al. 2019; Kremer 
et al. 2020b; Ye & Fragione 2022), and gravitational-wave sources 
made up of compact objects (Clausen et al. 2013; Rodriguez et al. 
2015, 2016; Askar et al. 2017; Arca Sedda 2020; Ye et al. 2020; 
Kremer et al. 2021). 


Depend on their distances from the Galactic Center (GC), the 
dynamical evolution of GICs is also subjected to the gravita- 
tional tidal field of the Milky Way (Baumgardt & Makino 2003). 
The expanding cluster halo could be truncated by the external 
tidal field, and tidal stripping may lead to the loss of stars from 
GICs, enhance the energy outflow of GICs and thereby 
accelerating the dynamical evolution of the cluster (Gnedin & 
Ostriker 1997; Gnedin et al. 1999). Besides, dynamical friction 
between GICs and Galactic background stars may lead to the 
spiral in of GICs toward the deep potential well of the Galaxy 
(Tremaine et al. 1975; Moreno et al. 2022), where tidal stripping 
would eventually lead to the fully dissolution of the clusters. 
GICs therefore can take their binary burning products into the 
Inner Galaxy (Fragione et al. 2018a; Arca-Sedda et al. 2018b), 
and contribute the cluster mass to the growing Galactic nuclei 
(Antonini et al. 2012; Antonini 2013; Gnedin et al. 2014). 

Thanks to their exclusive formation process and emission 
properties, the binary burning products are proved to be superb 
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tracing particles of GIC stellar dynamical interactions and 
cluster evolution. For example, compared with normal stars, the 
emission of LMXBs, MSPs, CVs, ABs and BSS are found to 
be much more luminous either in optical, X-ray, y-ray or radio 
band and can be detected and picked out from the dense core of 
GICs, making them ideal tracers for studying stellar dynamical 
encounters (Verbunt & Hut 1987; Pooley et al. 2003; 
Verbunt 2003; Pooley & Hut 2006), two-body relaxation and 
mass segregation effect of GICs (Ferraro et al. 2012; Cheng 
et al. 2019b, 20192). The cumulative cluster X-ray luminosity 
Lx and emissivity ex (Lx per unit stellar mass), proxies of 
population and abundance of weak X-ray sources (mainly CVs 
and ABs) in GICs (Cheng et al. 2018; Heinke et al. 2020), are 
also shown to be effective indicators for diagnosing cluster 
binary burning process (Cheng et al. 2018, 2020). 

Among all the binary burning products in GICs, MSPs are 
outstanding for many reasons. First, compared with BSS, CVs 
and ABs that can be formed through both primordial binary 
evolution and dynamical formation channels, it is widely 
recognized that LMXBs, the progenitors of MSPs, are 
dynamically formed in GICs (Verbunt & Hut 1987; Pooley 
et al. 2003; Verbunt 2003), and the abundances (number per 
unit stellar mass) of LMXBs and MSPs are more than 
7.100 times higher in GICs than in the Galactic field (Clark 
1975; Katz 1975; Camilo & Rasio 2005; Ransom 2008). 
Besides, unlike other tracers that could be contaminated by 
fore- and background sources, MSPs are mainly detected via 
either radio or y-ray observations, in which contamination is 
negligible. Moreover, the lifetime of MSPs is very long 
(~10'° yr) and their luminosities are observed to be very stable. 
Taking these aspects together, it is safe to say that MSPs are the 
best probe of cluster dynamical evolution, especially in tracing 
the tidal dissolution of GICs in the Milky Way. 

In y-ray astronomy, both MSPs and GICs are recognized as 
important y-ray emitters in the Galaxy (Abdo et al. 2010; Tam 
et al. 2011), and ~30 GICs have been identified as point sources 
in the fourth Fermi Large Area Telescope catalog (4FGL, 
Abdollahi et al. 2020). The y-ray emission of GICs is widely 
assumed to originate from MSPs that reside in the clusters 
(Cheng et al. 2010; Bednarek & Sobczak 2013), and the 
detection of pulsed gamma-ray emission from two GICs has 
further strengthened the case for this connection (Freire et al. 
2011; Wu et al. 2013). In concert with LMXBs, a strong 
correlation between ^-ray luminosity L., and the stellar encounter 
rate Г has been established among GICs (Abdo et al. 2010; Hui 
et al. 2011; Hooper & Linden 2016; Tam et al. 2016; Zhang 
et al. 2016; Lloyd et al. 2018; de Menezes et al. 2019), which 
lend a support to the dynamical origin of MSPs in GICs. 

On the other hand, the Fermi Large Area Telescope (Fermi- 
LAT) has discovered an unexpected y-ray excess around the GC, 
peaking at a few GeV, with an approximately spherical density 
profile xr 2“ and extending to 10°-20° (1.5-3 крс) from the GC 
(Goodenough & Hooper 2009; Hooper & Goodenough 201 1a; 
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Hooper & Linden 2011b; Di Mauro 2021). The nature of this 
"Galactic Center Excess" (GCE) is still not clear, and many 
promising models have been proposed over the past decade, 
including dark matter annihilation (Abazajian & Kaplinghat 2012; 
Daylan et al. 2016), emission of thousands of unresolved MSPs 
(Abazajian 2011; Yuan & Zhang 2014; Bartels et al. 2016), or 
emission from cosmic rays injected at ће GC (Carlson & 
Profumo 2014; Petrovié et al. 2014; Cholis et al. 2015; Gaggero 
et al. 2015). In the MSP scenario, Galactic GICs fallen into the 
GC are expected to deposit their population of MSPs, which are 
then inherited by the Galactic nuclear star cluster and nuclear 
bulge, and contribute to the observed y-ray excess (Bednarek & 
Sobczak 2013; Brandt & Kocsis 2015; Fragione et al. 20182; 
Arca-Sedda et al. 2018b; Abbate et al. 2018). 

In the present work, we perform a ^-ray study of the Galactic 
GICs based on the archival Fermi data. Unlike previous works 
that only focus on the cumulative GIC y-ray luminosity L, 
(Abdo et al. 2010; Hooper & Linden 2016; de Menezes et al. 
2019), we also measure the ^-ray emissivity €, (L, per unit 
stellar mass) of the GICs, and explore their dependence on 
various cluster parameters. The cumulative cluster ^-ray 
luminosity offers us a chance to study the population of MSPs 
at the Inner Galaxy, where single MSP is hard to be detected 
via radio or ^-ray observation. On the other hand, as 
demonstrated by Cheng et al. (2018) and Heinke et al. (2020) 
in the X-ray band, the emissivity has been proved to be a 
reliable indicator of exotic objects abundance in GICs, which is 
insensitive to the luminosity function (LF) and can be applied 
to a large cluster sample in a highly uniform fashion. More 
importantly, the derived GIC ^-ray emissivity can be directly 
compared to the stellar ^-ray emissivity of the Galactic nuclear 
star cluster and nuclear bulge, which are crucial to evaluating 
the relative abundance of putative MSPs in these environments, 
thus estimating the possible stellar origins and contribution to 
the GCE by dissolved GICs. 

The limitation of our approach is that we assume the 
measured cluster y-ray luminosity is a good proxy of the 
population of hosted MSPs, that may not be true for all clusters, 
especially for GICs contains few MSPs. However, as we will 
show below, although the uncertainty of small counts of MSPs 
in individual GICs may introduce large scatter to our cluster 
sample,” but are unlikely to create the correlations and trends 
observed in this work. 

The paper is organized as follows. Section 2 describes the 
data reduction and analysis that lead to the detection and the 
measurement of the y-ray luminosity L, and emissivity є; of 
individual GICs. Section 3 explores the dependence of L, and 
€, on various cluster physical properties. A discussion and a 
summary of our results are presented in Sections 4 and 5, 


5 Ih fact, the scatter of derived correlations in this paper is about an order of 
magnitude larger. 
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respectively. Throughout this work we quote 1c errors, unless 
otherwise stated. 


2. -y- Ray Data Analysis 
2.1. Data Reduction and Analysis 


We analyzed the archival Fermi-LAT data of the 157 GICs 
presented in the catalog of Harris (1996). The Fermi data was 
observed from 2008 August 8 to 2020 November 8 (MET: 
239 846 401-626486405), with a time span of ~12 yr. We 
used the Fermi tools release 2.0 to analyze the data, with the 
energy band was restricted to 100MeV—300 GeV and divided 
into 15 logarithmically spaced energy bins. We selected the 
events with source class (evclass — 128, evtype — 3) and filter 
the data with DATAQUAL > 0, LATCONFIG == 1. A zenith 
angle cut of <90° and a satellite rocking angle cut of <52° was 
applied to avoid contamination from the Earth limb. 

The region of interest (ROI) of each target was restricted to a 
14? x 14? rectangular box centered on the optical center of the 
GICs. We used the Make4FGLxml.py tool and the 4FGL DR2 
catalog to generate the background source list within the ROI. For 
diffuse background modeling, we adopted the most recent Galactic 
interstellar emission model gll iem v07.fits and the isotropic 
spectral template iso PSR3 SOURCE V3 vl1.txt. The instrument 
response function (IRF) was set as PSR3. SOURCE УЗ. 

All targets were investigated by means of binned likelihood 
analysis (gtlike tooI—DRMNFB, NEWMENUIT algorithm). 
To search for the y-ray emission from the GICs, we added а 
putative point source at the optical center of the cluster, and 
used the gttsmap tool to derive the TS maps. From the TS maps, 
one can visually inspect the detections of unkonwn sources 
located within the ROI. We then removed the new sources by 
adding them to the background source lists and refit the data 
again. The spectral models of GCs were set as either a power- 
law (PL), a PL plus exponential cut-off (PL--Expcut), or a 
logparabola (LP) model. During the fitting, all the GIC 
parameters (i.e., coordinates, photon index, cut off energy 
and normlization) were set as free, whereas for diffuse 
backgrounds and point sources located within 5? of the ROI, 
only the normalization parameter was left free to vary. To 
quantify the significance of the detection, a test statistic (TS) 
was calculated, defined as TS = 2(log Lı — log Lo), where Г 
(Lo) is the maximum likelihood value of the model with 
(without) the putative source. The chosen criteria for detection 
was TS > 25, corresponding to a significance slightly above 4o. 


2.2. Data Analysis Results 


For the 157 GICs tabulated in the catalog of Harris (1996), 
about 2596 (39/157) of them are found to be 4-ray bright 
(TS > 25) in this work, with two clusters, HP 1 (TS = 44.3) and 
Terzan 9 (TS = 42.1), are first identified as y-ray emitters. Cross- 
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checking the 39 clusters with the 4FGL-DR3 catalog (Abdollahi 
et al. 2022), five clusters, NGC 362 (TS = 16) and NGC 6304 
(TS = 21.7), are found to have a smaller significance than our 
detection threshold. Besides, Liller 1 was found to be y-ray 
bright by Tam et al. (2011), while our data fitting suggests that 
the detection is marginal (TS = 21.6). All the 39 GICs are 
located within a distance of D=15kpc from the Earth 
(Figure 1(a)), and there is no evident observation bias for GICs 
near the GC (Figure 1(b)) and the Galactic disk (Figure 1(c)), or 
significant dependence on cluster metallicity (Figure 1(d)). 
However, it seems that GICs are more likely to be identified as 
y-ray emitters, provided that they are more massive (Figure 1(e)) 
or have larger stellar encounter rate (Figure 1(f)). Interestingly, 
the y-ray detection rate of core-collapsed GICs (12/29) is found 
to be ~2 times higher than the dynamically normal GICs (27/ 
128), which suggest that the dynamical evolution history is also 
an important factor influence the y-ray luminosity of GICs. 

The y-ray data analysis results are summarized in Table 1, 
which is segmented into two panels according to the core 
collapse flag in the catalog of Harris (1996). The likelihood fit 
parameters arranged from left to right are: cluster name, 
spectral models, test statistic values, PL photon index, 
exponential cut-off energy E. and measured GIC energy flux 
fy Adopting the cluster distance (D) presented in the catalog of 
Baumgardt & Hilker (2018), we also calculated the cumulative 
GIC ^-ray luminosity with function L, = 4nDf.. We define the 
GIC ^-ray emissivity as €, — L,/M, with M being the cluster 
mass given in the catalog of Baumgardt & Hilker (2018). L, 
and є, are listed in Columns (7) and (8) of Table 1. The errors 
are given in lo standard deviation level. 

To check the spatial consistency between the ^-ray emission 
and the optical centers of GICs, we plot in Figure 2 the TS maps 
of the 39 clusters. For each GIC, the image was restricted to a 
5? x 5? box, with the central green circle indicates the tidal 
radius of the cluster. All maps have evidence for good 
coincidence between the ^-ray emission and the GIC centroids, 
except for NGC 1904, where the offset between the optical 
center and the peak of the y-ray emission (cayan ellipse) is 
~0°3. Besides the coincidence of spatial coordinates, it is 
possible that the y-ray emission is associated with a fore/ 
background source rather than the cluster. We estimated this 
probability with function Pandom = aR2N /Sgor, where R, is 
the GIC tidal radius, Spo; is the area of the ROI, and N is the 
number of ^-ray sources detected within Spo;. The values of 
Pandom аге presented in the last column of Table 1. In most 
clusters, the Pandom value is less than ~10%, the exceptions are 
NGC 104 (16.9%), NGC 5139 (34.5%), NGC 6624 (18.1%), 
NGC 6752 (32.1%) and NGC 6656 (40.1%), where the random 
detection probability is over ten percent However, as 
illustrated in Figure 2, А, of these clusters are also found to 
be very large, thus the larger Prandom is more likely to result 
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Figure 1. Histogram distributions of GICs as a function of cluster parameters: distance from the Earth (a), distance from the Galactic center (b), distance from the 
Galactic disk (c), cluster metallicity (d), cluster mass (е) and the stellar encounter rate (f). The blue lines represent the y-ray bright GICs, while clusters without ^-ray 


detection are shown as green lines. 


from the larger tidal radius rather than being associated with a 
fore/background source. 


3. Statistic Relations 


As discussed in Section 1, the measured GIC 4-ray 
luminosity (L5) arises from a collection of MSPs reside in the 
cluster. Generally speaking, there are two dominant factors 
may influence the number of MSPs of a given GIC: (1) the 
population (thus the cluster mass M) of neutron stars (NS) 
hosted by the cluster, and (ii) the stellar dynamical interactions 
that transfrom NS into LMXBs and MSPs. The second factor 
also refers as the binary burning process, which is thought to be 
the essential internal energy source of cluster dynamical 
evolution. In this regard, the y-ray emissivity (є, = L.,/M) is 
better than L, in tracing the formation efficiency of MSPs in 
GICs, and both L, and e, may serve as sensitive probe to study 
the dynamical evolution of GICs. With this in mind, we 
examine the dependence of L, and c, on various cluster 
physical parameters in this section. If tidal stripping have 
played an important effect in the dynamical evolution of GICs, 


significant correlations between GIC ^-ray emission and the 
Galactic environment parameters are also expected. The GIC 
parameters are mainly taken from Baumgardt & Hilker (2018) 
or Harris (1996), unless the specific references are quoted. 
Throughout the paper, we also classified the GICs into two 
subgroups (i.e., the dynamically normal and core-collapsed 
GICs) according to their labels in Harris (1996). The core- 
collapsed GICs are generally considered to be in the late phase 
of cluster dynamical evolution, and thus are dynamically more 
older than the normal GICs. 


3.1. Correlations with Cluster Mass 


In Figure 3 (a), we first explore L, versus M for all GICs. The 
GIC mass is adopted from the online catalog? of Baumgardt & 
Hilker (2018), which is derived by comparing the observed 
cluster velocity dispersion, surface density, and stellar mass 
function profiles against N-body simulations. Although the 
scatter is substantial, a moderate positive correlation between L, 


6 https: / [people.smp.uq.edu.au/HolgerBaumgardt/globular/ 
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and M is observable, from which we find the Spearman's rank 
coefficient r — 0.454 and r — 0.464 for the total and dynamically 
normal GICs, with p = 0.004 and p = 0.017 for random 
correlation. We fit a power-law function to the correlations and 
obtained Lg and Іо M^? 14029 for total and the 
dynamically normal GICs, respectively. The best fitting func- 
tions are shown as the blue and the olive lines in Figure 3 (a), 
and the dotted curves mark the 95% confidence range. 
Compared with the dynamically normal GICs, the dependence of 
L., on M is marginal for core-collapsed GICs, with Spearman's 
rank coefficient r — 0.259 and random correlation p-value of 
р = 0.417. The minimum and median y-ray luminosity of 
dynamically normal GICs is measured to be Ly min = 5.41 x 
1023 erg s~! and Lymea=4.1 x 10°* erg s !, which is compar- 
able to that of the core-collapsed GICs (i.e., Ly min = 3.45 x 
103% erg s~! and Г. „оа = 3.22 x 10° ergs !). The exception is 
the most luminous GICs (ie, L, Z 10% ergs '), they are 
dominated by dynamically normal GICs and thus are responsible 
for the much larger scatter of L, in dynamically normal GICs 
than the core-collapsed ones. 

In Figure 3 (b), we examine the dependence of €, on M. It can 
be seen that c, has a substantial scatter ranging from 
~10% еге 5! M;! to ~10% erg s! M;!. Nevertheless, a 
marginally significant negative correlation between є; and M is 
suggested by the Spearman's rank coefficient, with r — — 0.396, 
r= — 0.343 апа r= — 0.359, and random correlation p-value 
of p — 0.045, p — 0.276 and p — 0.027 for the dynamically 
normal, core-collapsed and total sample of GICS, respectively. A 
power-law function fitting gives €, x M- 0495022 (olive lines) 
and cam ee (blue lines) for the dynamically normal 
and total GICs. These negative correlations are consistent 
with the sub-linear correlations found in Figure 3 (a) The 
average ^-ray emissivity and standard deviation of all GICs is 
measured to be (є,) = (2.48 + 0.41) x 10? erg s^! M5' and 
ог = 2.52 x 10? erg s^! Мс', and no significant difference 
between the dynamically normal GICs ((&) = (2.58 + 0.54) х 
10? erg s"! M5!, о, = 2.75 x 10? erg s^! М!) and core- 
collapsed GICs ((€,) = (2.25 + 0.58) x 10? erg s^! M5, о, = 
2.02 x 10? erg s^! М!) are found. 

The dependence of L, on M is consistent with the work of 
Hooper & Linden (2016), where massive GICs are also found 
to have larger y-ray luminosities. This trend is also in 
agreement with Figure 1 (е), supporting a y-ray detection 
preference for massive GICs. However, the negative є, — M 
correlation suggests that the cluster mass is not the decisive 
factor influencing the abundance (or formation efficiency) of 
MSPs in GICs. We note that a similar e, — M relationship has 
also been reported by Fragione et al. (20182) in their Figure 3. 
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3.2. Correlations with Stellar Encounter Rates 


Traditionally, there are two parameters in literature to 
describe the stellar encounter rate of GICs. The first is the 
total stellar encounter rate, Г c /о?/ c, an integral of stellar 
encounters over the cluster volume, with p the stellar density 
and c the cluster velocity dispersion (Verbunt & Hut 1987; 
Verbunt 2003). The second is the specific encounter rate, 
defined as Л = Г/М, and measure the chance that a particular 
star (or binary) undergoes dynamical encounters in GICs 
(Pooley & Hut 2006). For this reason, Г is thought to be a 
resonable indicator of the total amount of exotic objects that 
can be dynamically produced in GICs (Pooley et al. 2003; 
Cheng et al. 2018). While A is considered to be highly 
correlated with the abundance (or the formation efficiency) of 
exotic objects in GICs (Pooley & Hut 2006; Cheng et al. 2018). 

In Figure 4 (a), we plot L, versus Г for each GIC. The value 
of Г are adopted from Bahramian et al. (2013), which is 
integrated over the cluster and normalized to a value of 1000 
for NGC 104 (Table 2). It is clear that L, is highly correlated 
with Г for dynamically normal GlCs. The Spearman’s 
correlation coefficient is r = 0.697, with a random correlation 
p-value of p= 1.1 x 107“. If core-collapsed GICs are taken 
into account, this correlation is still significant, with r — 0.583 
and р = 1.5 х 10 ^ for total GICs. The best-fitting functions 
can be written as «Г олко (olive curves) and 
Lyx [973:0.!! (blue curves) for the dynamically normal and 
total GICs, respectively. These correlations are consistent with 
the finding of de Menezes et al. (2019), where the possible 
number of MSPs (N) of GICs is found to be proportional to Г, 
with N x D?9*59-15 in a sample of 23 GICs. 

According to the Spearman’s rank correlation coefficients, it is 
clear that L, depends on I' better than M, which suggest that 
stellar dynamical interaction is the fundamental factor that 
influence the population of MSPs in GICs. To minimize the 
dependence on cluster mass, we also examine c, versus A in 
Figure 4 (b) The value of A was calculated with equation 
А —T/Me, with Mg be the cluster mass in units of 106M.... Since 
I' was estimated based on the V-band luminosity (Bahramian 
et al. 2013), here we calculated Mg using the V-band-based 
magnitude for consistency, following an empirical mass-to- 
magnitude relation described in Cheng et al. (2018). From 
Figure 4 (b), it is obvious that the dynamically normal GICs have 
a larger €, with increasing A. The Spearman's rank correlation 
coefficient is r = 0.616, with p=1.0 x 10? for random 
correlation. This correlation becomes less significant when taking 
the core-collapsed GICs into acount, with r = 0.368 and 
p = 0.025 for the total GICs. A power-law fitting yields a function 
of e, x 978:9.5 (olive curves) for the dynamically normal GICs. 

Compared with dynamically normal GICs, the core- 
collapsed GICs are found to have much large scatter in 
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Table 1 
^-Ray Data Analysis Results of the 27 Dynamically Normal GICs (Upper Panel) and 12 Core-collapsed GICs (Lower Panel) Detected in this Work 
Name Model TS Photon Index E. ГА L, €, Prandom 
(1) (2) (3) (4) (5) (6) (7) (8) (9) 
Dynamically Normal GICs 
NGC 104 PL+Expcut 5557.3 1.39 + 0.07 2.82 + 0.30 26.9 + 0.7 65.9 + 1.8 7.37 £0.21 16.9% 
NGC 1851 PL+Expcut 42.8 1.46 + 0.24 1.71 + 0.32 1.25 + 0.32 18.2 + 4.9 5.71 + 1:53 0.596 
NGC 2808 PL+Expcut 102.5 1.22 + 0.37 2.29 + 1.00 2.99 + 0.47 36.3 + 5.7 4.20 + 0.66 0.896 
NGC 5139 PL+Expcut 977.3 1.20 + 0.17 1.98 + 0.37 11.6 + 0.8 41.0 + 2.8 1.13 + 0.13 34.596 
NGC 5286 PL 25.4 2.67 + 0.19 ө 3.11 + 0.67 46.0 + 9.9 13.0 + 2.8 0.8% 
NGC 5904 PL+Expcut 46.1 0.93 + 0.7 2.22 + 1.30 1.18 + 0.26 7.92 £1.77 2.01 + 0.45 4.8% 
NGC 6093 PL+Expcut 105.6 1.94 + 0.2 6.29 + 2.10 3.84 + 0.72 49.3 + 9.3 14.6 + 2.8 0.596 
NGC 6139 PL+Expcut 63.5 1.94 + 0.25 4.09 + 2.24 5.05 + 0.90 61.1 + 10.8 18.9 + 3.9 3.8% 
NGC 6218 PL+Expcut 30.3 3.5+0.01 0.10 + 0.00 2.54 + 0.58 7.96 + 1.82 7.44 + 1.72 2.5% 
NGC 6254 LP 57.4 © ө 1.98 + 0.30 5.43 + 1.01 2.65 + 0.50 2.996 
NGC 6316 PL+Expcut 309.0 1.93 + 0.14 5.32 + 1.87 12.8 + 1.1 190.9 + 16.3 60.0 + 9.0 3.096 
NGC 6341 PL+Expcut 59.9 1.89 + 0.62 2.93 + 2.56 2.48 + 0.36 21.5 + 3.1 6.11 + 0.88 2.396 
NGC 6388 PL+Expcut 1266.5 1.46 + 0.11 2.45 + 0.39 19.2 + 0.4 287.4 + 6.1 23.0 + 0.5 1.9% 
NGC 6402 PL+Expcut 56.0 1.51 + 0.38 3.71 + 1.89 2.80 + 0.62 28.1 + 6.2 4.74 + 1.07 0.7% 
NGC 6440 PL+Expcut 321.6 2.09 + 0.12 8.78 + 3.29 20.3 + 1.3 165.8 + 10.5 33:9: 3.0 1.796 
NGC 6441 PL+Expcut 551.1 1.71 + 0.09 3.03 + 0.34 21.6 + 1.1 420.0 + 21.0 31.8 + 1.6 3.096 
NGC 6528 PL+Expcut 40.9 1.29 + 0.26 1.55 + 0.31 3.52 + 0.89 25.9 + 6.6 45.7 + 12.7 1.1% 
NGC 6626 PL+Expcut 1102.2 1.51 + 0.05 1.16 + 0.06 ААЛ ЕТ 73.0 + 3.8 24.4 + 2.0 6.296 
NGC 6637 LP 29.1 © es 1.49 + 0.35 14.2 € 3.3 9.14 + 2.15 2.396 
NGC 6652 PL+Expcut 158.7 1.49 + 0.4 2.21 + 1.39 4.84 + 0.60 52.0 + 6.4 108 + 20.6 1.2% 
NGC 6656 PL+Expcut 100.5 0.83 + 0.29 1.05 + 0.13 4.17 + 0.47 5.41 + 0.85 1.14 + 0.18 40.196 
NGC 6717 PL+Expcut 86.1 0.63 + 0.52 1.82 + 0.70 2.50 + 0.43 17.0 € 2.9 47.4 + 13.6 1.896 
NGC 6838 PL 31.8 2.71 + 0.2 et 3.92 + 0.79 7.53 + 1.51 16.5 + 3.4 1.5% 
2MASS-GCOI PL+Expcut 92.8 1.87 + 0.19 3.55 + 1.29 19.6 + 2.1 26.7 + 2.9 76.3 + 8.3 1.8% 
Glimpse 01 PL+Expcut 760.2 1.65 + 0.16 4.13 + 0.83 78.2 + 2.7 108.5 + 3.8 36.2 + 3.8 6.3% 
Glimpse 02 PL+Expcut 138.1 1.47 + 0.13 1.62 + 0.18 30.1 + 2.4 76.4 + 6.0 et 10.196 
Terzan 5 PL+Expcut 5290.6 1.74 + 0.03 3.98 + 0.18 125.0 + 2.2 657.3 + 11.5 70.3 + 5:3 2.996 
Core-Collapsed GICs 

NGC 1904 PL 53.9 2.54 + 0.16 ө 2.5 + 0.4 52.3 + 8.5 37.7 + 6.8 0.6% 
NGC 6266 PL+Expcut 1332.0 1.45 + 0.11 2.75 + 0.44 19.4 + 0.8 95.6 + 3.9 15.7 + 0.7 5.2% 
NGC 6397 PL+Expcut 44.8 2.15 +055 3.89 + 4.09 4.67 + 0.68 3.45 + 0.50 3.57 + 0.52 6.4% 
NGC 6541 PL+Expcut 136.1 1.61 + 0.34 2.89 + 1.65 4.51 + 0.58 31.3 +4.1 10.7 + 1.4 4.1% 
NGC 6624 PL+Expcut 812.1 1.15 + 0.22 1.22 + 0.44 13.10 € 1.1 101.1 + 82 64.8 + 5.5 18.196 
NGC 6723 PL+Expcut 29.8 1.54 + 0.43 4.71 + 5.73 1.75 + 0.48 14.4 + 3.9 8.11 + 228 1.796 
NGC 6752 PL+Expcut 193.1 0.59 + 0.43 1.06 + 0.28 3.26 + 0.39 6.64 + 0.79 2.41 + 0.29 32.1% 
NGC 7078 PL+Expcut 103.0 1.98 + 0.4 1.55 + 0.92 3.64 + 0.53 50.1 + 7.2 7.9] + 1.14 7.096 
HP1 PL+Expcut 44.3 0.48 + 0.24 1.55 + 021 5.62 + 0.98 33.0 + 5.8 26.6 + 5.9 5.696 
Terzan 1 LP 78.0 Ut © 6.34 + 1.05 24.5 + 4.1 16.3 + 3.8 10.6% 
Terzan 2 PL+Expcut 58.4 0.55 + 0.22 1.55 + 0.2 10.4 + 1.5 75.0 + 10.4 55.1 :E 12.7 5.596 
Terzan 9 PL+Expcut 42.1 1.47 + 0.35 1.55 + 0.26 6.28 + 3.94 25.1 + 15.7 20.9 + 13.3 6.0% 


Note. Columns 1-2: Name of GICs and Spectral models. Columns 3—4: The test statistic value and the power-law phonton index. Columns 5—7: Cut-off energy (in 


units of GeV), GIC y-ray flux (in units of 10 ^ erg ст 257! 


1028 erg s-! M ^! C) and the random detection probability. 


Figure 4. The difference may arises from the rapid change of 
cluster structure in these systems. Once GICs are running out of 
binary systems and undergo deep core collapse, their structure 
are not stable any more and the cluster core may suffer from 
gravothermal oscillations (Fregeau et al. 2003), core-collapsed 
GICs therefore may create extremely frequent stellar encounter 
environments. Indeed, some core-collapsed GICs (e.g., NGC 


) and luminosity (in units of 10? erg s). Columns 8-9: GIC y-ray emissivity (in units of 


6624 and NGC 7078) are proved to be clusters with the highest 
stellar encounter rates in Table 2, supporting a currently deep 
core collapse in these systems. While some clusters (e.g., HP 
1 and Terzan 1) are found to have the smallest values of Г 
and A, which may indicate a core bounce phase after the deep 
core collapse. On the other hand, with extremely large A, it is 
possible for LMXBs to be dynamically disrupted during deep 
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Figure 2. (a) 5? x 5? test statistic maps of the 39 GICs with TS > 25 in Table 1. The GIC name is shown as green text in the upper right of each panel, and the green 
circle represents the tidal radius of the cluster. The y-ray emission are in coincident with the optical centers of the GICs. The exception is NGC 1904, where the la 
error (magenta circle) of the y-ray emission peak is in coincident with the 95% error ellipse (cyan) of the 4FGL catalog, which has an offset of ~0°3 from the cluster 
center. (b) Continued test statistic maps of the 39 GICs presented in Table 1. 
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Figure 3. GIC ^-ray luminosity (a) and y-ray emissivity (b) as a function of cluster mass. The olive and red pluses represent the dynamically normal and core- 
collapsed GICs separately. The red, olive and blue texts donate the Spearman’s rank correlation coefficient of the core-collapsed, dynamically normal and the total GIC 
samples, respectively. The solid lines and associated dotted curves are the best fitting functions and the 95% confidence range: red for core-collapsed GICs, olive for 


dynamically normal GICs, and blue for the total sample. 


core collapse (Verbunt & Freire 2014), and strong encounters 
may lead to the ejection of MSPs from the clusters 
(Section 4.2). These processes may introduce uncertainties 
to the L,—I and e,— А correlations for core-col- 
lapsed GICs. 


3.3. Correlations with Cluster Structure Parameters 


To find out the dependence of GIC y-ray esmission on GIC 
dynamical evolution history, we sort the cluster structure 
parameters and investigate their correlations to L, and є, in 
Figures 5 and 6. Generally speaking, MSPs are prone to be 
dynamically formed in the dense core of GICs. Therefore, we 
first plot L, versus GIC core radius R., core density pe, core 
relaxation time f,, and central velocity dispersion с, in the upper 
panels of Figure 5. No clear correlation exists for the total GIC 
sample, except that L, is highly correlated with о, (r = 0.574, 
р= 2.0 х 10 5) in Figure 5(d). When only the dynamically 
normal GICs are considered, this correlation is still evident, with 
Spearman's r — 0.672 and random correlation p-value of 
p —24 x 10 ^. We fit this correlation with power-law function, 
which gives L, x g239*9^^ (blue lines) and Ly x аав 
(olive lines) for the total and dynamically normal GICs, 
respectively. However, it seems that the cluster mass is the 
more fundamental parameter underlying the Г.с, relation, since 
the dependence of €, on о, is less significant in Figure 6(d), and 
с. is highly correlated with M in GICs (Table 3). 

Unlike the linear relations reported Sections 3.1 and 3.2, the 
Spearman's rank correlation coefficients suggest that the 


dynamically normal GICs and core-collapsed GICs аге 
statistically different in Figures 5(a)-(c). As GICs evolve to 
older dynamical age (1.е., with decreasing /,„) and become more 
compact (i.e., with smaller R, and larger p,), the dynamically 
normal GICs are more likely to exhibit larger L,. Whereas for 
core-collapsed GICs, this trend is reversed. When the influence 
of cluster mass is eliminated, this tendency still holds in 
Figures 6(a)-(c), where the dynamically normal GICs are 
measured to have larger є, with increasing dynamical age, 
which then reverses to smaller values as clusters evolve into 
core-collapsed GICs. We suggest that these features may 
indicate an episodic ejection of MSPs from GICs. Namely, 
during the deep core collapse and subsequent core bounce, 
significant number of MSPs may have been dynamically 
ejected from the host clusters. The implication of this effect 
will be addressed in Section 4.2. Here, we find the significant 
correlations of dynamically normal GICs can be expressed as 
L, x о e, x R; 0955027 and €, о 10613 

plotted as olive lines in the figures. 

In the middle panels of Figures 5 and 6, we test the 
dependence of L, and є; on the GIC half-mass-radius Ку, 
stellar density (оу) inside А, half-mass relaxation time fp and 
central escape velocity v... It is obvious that the dynamically 
normal GICs have a larger L, with increasing p, and vesc, and 
the best-fitting functions (olive lines) can be written as 
Ly x T апа L, x v215596! Tf core-collapsed GICs 


-0.11 They are 


esc 


are included, these correlations are still evident, with 
L 0.58+0.13 2.25+0.48 
yX p, and Ly c уг for the total GICs (blue 
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Figure 4. GIC y-ray luminosity as a function of the stellar encounter rate Г (a) and emissivity as a function of the specific encounter rate A (b). The color-coded 
symbols denate different types of GICs, as in Figure 3, and the same color coded text gives the Spearman's rank correlation coefficients: red for core-collapsed, olive 
for dynamically normal, and blue for the total. The olive and blue solid lines mark the best-fitting funtion of the dynamically normal and the total GICs, and dotted 


curves represent the 95% fitting confidence range. 


lines). However, the dependence of e, on А, and Vege аге not 
significant in Figures 6(f) and (h), reflecting again that cluster 
mass is the more fundamental parameter underlying these 
relations. On the other hand, GICs with smaller А, appears to 
have larger L,, as indicated by the mild anti-correlation 
coefficients in the figures. This trend is still evident in 
Figure 6(e), where the best fitting functions between c, and 


К, can be expressed as e, с jg BL Ey Ren and 


€, X Ry 1765048 for the dynamically normal, core-collapsed 
and total samples of GICs. Besides, although the dependence of 
L, on fj, is not apparent in Figure 5(g), €, is found to decrease 
with inceasing fn in Figure 6(g), with the best fitting functions 


can be written as e, сс po €, x fi n aad €, Xx 


150802017 for dynamically normal, core-collasped and total 


GICs, respectively. 

By and large, the dependence of y-ray esmission on GIC 
half-mass paramters is in agreement with that on the core 
paramters. That is, as GICs evolve to older dynamical age 
(tj x tp! or tg x th) and become more compact, GICs are 
more likely to exhibit larger L, and є,. The discrepancy may 
arise from the core collapse, in which the cluster core structure 
may change dramatically during the gravothermal oscillations. 
Whereas for the half-mass parameters, their response to core 
collapse will be more moderate (Fregeau et al. 2003), thus both 
dynamically normal and core-collapsed GICs are found to 
exhibit similar correlations in Figures 5(e)-(h) and 6(e)-(h). 

In the bottom panels of Figures 5 and 6, we examine the 
dependence of L, and ce, on the cluster tidal radius R, stellar 


10 


mass function (MF) slope o, concentration parameter c and 
metallicity [Fe/H]. L, is found to be independent of R,, a and c 
in Figures 5(1), (j) and (k). However, It is clear that GICs with 
smaller R, and larger o are more likely to have larger c, 
(Figures 6(1) and (j). The power-law fitting yields 
€, о К, 1085028. €, x RyC098:037 ang e, x 190022 for 
the dynamically normal, core-collapsed and total GICs, 
respectively. While for the €,—a relation, which can be expressed 
as loge, = аа +b in Figure 6(j, with the best-fitting 
parameters а = (0.57 + 0.10) and b = (29.24 + 0.07) for total 
GICs, a = (0.67 + 0.14) and b = (29.26 + 0.09) for dynamically 
normal GICs, and a= (0.41 + 0.15) and b = (29.22 + 0.11) for 
core-collapsed GICs, respectively. 

Observationally, many authors argued that LMXBs are more 
likely to be detected in old, metal-rich GICs, rather than in 
young, metal-poor ones (Sivakoff et al. 2007; Li et al. 2010; 
Kim et al. 2013). As the offspring of LMXB, Figure 5(1) also 
suggests that MSPs are more likely to be found in metal-rich 
GICs. We fit the L,-[Fe/H] relation with the function 
log L, = a[Fe/H] + b, the  best-fitting parameters are 
а = (0.53 + 0.18) and b= (35.16 + 0.21) for dynamically nor- 
mal GICs, а= (0.47 + 0.14) and b = (35.09 + 0.17) for total 
GICs, respectively. When the influence of cluster mass was 
eliminated, we find this tendency are still evident in Figure 6(l). 
The best-fiting function log є; = a[Fe/H] + b gives para- 
meters of a= (0.58 + 0.17) and b = (29.71 + 0.20) for dynami- 
cally normal GICs, a = (0.54 + 0.18) and b = (29.88 + 0.26) for 
core-collapsed GICs, and а = (0.53 + 0.13) and b= (29.73 + 
0.13) for total GICs, respectively. We note that a similar 
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Table 2 
Integrated Stellar Ecounter Rate (Г) and Specific Encounter Rate (A) of GICs 
Name r —ó +6 A —ó +6 
Dynamically Normal GICs: 
NGC 104 1000 134 154 844 113 130 
NGC 1851 1530 186 198 3524 428 456 
NGC 2808 923 83 67 801 72 58 
NGC 5139 90.4 20.4 26.3 35.2 7.9 10.2 
NGC 5286 458 61 58 723 96 92 
NGC 5904 164 30 39 243 45 57 
NGC 6093 532 69 59 1344 174 149 
NGC 6139 307 82 95 688 184 214 
NGC 6218 13.0 4.0 5.4 76.6 23.8 32.1 
NGC 6254 31.4 4.1 4.3 158 21 22 
NGC 6316 77.0 14.8 25.4 176 34 58 
МОС 6341 270 29 30 695 75 TI 
NGC 6388 899 213 238 766 181 203 
NGC 6402 124 30 32 141 34 36 
NGC 6440 1400 477 628 2190 746 983 
NGC 6441 2300 635 974 1600 442 678 
NGC 6528 278 50 114 3239 577 1328 
МОС 6626 648 91 84 1746 245 226 
МОС 6652 700 189 292 7508 2027 3132 
МОС 6656 77.5 25.9 31.9 153 51 63 
NGC 6717 39.8 13.7 21.8 1072 369 587 
NGC 6838 1.47 0.14 0.15 41.5 4.0 4.2 
2MASS-GCO01 n E us iis n m 
Glimpse 01 400 200 200 8560 4280 4280 
Glimpse 02 ET 25 S oe m sis 
Terzan 5 6800 3020 1040 3400 1510 520 


Core-Collapsed GICs: 


NGC 1904 116 45 68 412 159 240 


NGC 6266 1670 569 709 1758 599 741 
NGC 6397 84.1 18.3 18.3 919 200 200 
NGC 6541 386 63 95 746 122 184 
NGC 6624 1150 178 113 5742 889 564 
NGC 6723 11.4 4.3 8.0 41.6 16.0 292 
NGC 6752 401 126 182 1605 504 729 
NGC 7078 4510 986 1360 4705 1029 1419 
НРІ 0.66 0.30 0.41 8.51 3.87 5.29 
Terzan 1 0.29 0.17 0.27 24.7 14.5 23.0 
Terzan 2 22.1 14.4 28.6 486 317 629 


Terzan 9 1.71 0.96 1.67 278 156 271 


Note. The lower and upper limits are given at lo level. 


L,-[Fe/H] relation was also obtained by de Menezes et al. 
(2019) with 23 GICs. 


3.4. Correlations with the Galactic Environment 


Contrary to the contraction of the central regions, the 
outskirts of GICs are driven to expand to be truncated by the 
tidal force of the host galaxy. The tidal radius, defined as the 
boundary where a GIC may lose its stars to the Galaxy, is 
subject to its coordinates in the Milky Way (Baumgardt & 
Makino 2003) R, = (GM/2W3) R. Here С is the 


Feng et al. 


gravitation constant, Vg the circular velocity of the Galaxy and 
К.с the distance of the cluster from the GC. The independence 
of L, on R, in Figure 5(1), and the negative c,-R, relation in 
Figure 6(j), suggests that tidal stripping effect is another 
important factor influence the abundance of MSPs in GICs. 
Indeed, this tendency also can be confirmed by the positive 
eya correlation, since mass segregation and tidal stripping 
effects have a preference for the depletion of low-mass stars 
from GICs, thus tidally stripped clusters are more likely to 
exhibit larger stellar MF slope a (Vesperini & Heggie 1997; 
Baumgardt & Makino 2003). 

Therefore, we also investigate the dependence of GIC ^-ray 
emission on the Galactic environment parameters in Figures 7 
and 8. The coordinates of GICs are adopted from the online 
catalog of Baumgardt et al. (2019), with X point from the 
Galactic center toward the Sun, Y point in the direction of 
Galactic rotation at the Solar position, and Z point toward the 
North Galactic pole, respectively. In Figures 7(a) and 8(a), both 
L, and e, are found to decrease moderately with increasing X. 
However, these features may reflect an observational selection 
effect of GIC ^-ray detection by Fermi-LAT at the position of 
Sun (i.e., with X z 8 kpc, Y= 0 kpc and Z — 0 kpc), rather than a 
physical dependence of GIC ^-ray emission on X. Nevertheless, 
it seems that L, tend to reach the maximum values near the GC 
in Y and Z direction (Figures 7(b) and (c)), and this tendency is 
still evident when the influence of cluster mass was considered 
in Figures 8(b) and (c). We further plot in Figures 7(d)-(f) L, as 
a function of cluster distance to the GC (К), the perigalactic 
distance (К.с рег) and apogalactic distance (А, аро) of GIC orbit 
in the Milky Way. Although with large scatter, the mild anti- 
correlations between L, and the GIC distances to GC is in 
agreement with the finding of Figures 7(b) and (c), suggesting 
that more MSPs tend to be dynamically formed in GICs closer to 
the GC. This tendency is even more significant in Figures 8(d)- 
(f) where the best-fitting functions can be expressed as 


—1.13+0.21 —0.74£0.13 —1.21+0.20 
€, X Ry s Ey X Roo per and e, сс Ree apo for total 
—143::028 —0.73+0.16 —1.524027 
GICs, e, ос К ‚ Ey X Ryo per and e, сс Ree apo 


for dynamically normal GICs, and e, ос 07034 


RI aaa and є; х т m for core-collapsed GICs, 
respectively. 

For comparison, it will be interesting to mark in Figures 8(d)— 
(f) the 5-ray emissivity of stars in the Inner Galaxy, since tidally 
disrupted GICs are thought to have a contribution to at least part 
of the stars therein. With a total mass of (1.4 0.6) x 10? Ms 
for nuclear bulge stars and (0.91 + 0.7) x 10'°M 5 for boxy 
bulge stars, Bartels et al. (2018) showed that the emssion profile 
of GCE can be better fitted by stellar mass in the boxy and 
nuclear bulge rather than the dark matter profiles, and the y-ray 
emissivity of the Galactic boxy bulge and nuclear bulge is 
measured to be є; = (1.9 + 0.2) x I0? erg s^! М5! and 
€, = (1.1 + 0.6) x 10” erg s^! М! separately. We mark the 
boxy bulge and nuclear bulge with blue and olive rectangles in 


‚ &y © 
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Figure 5. Dependence of the y-ray luminosity on various physical properties of GICs. Top panels, from left to right: cluster core radius R,, core density log(p.), core 
relaxation time log(t;-), and central velocity dispersion o,. Middle panels, from left to right: cluster half-mass-radius R}, density inside half-mass-radius log(p,), half- 
mass relaxation time log(t,,), and central escape velocity vesc- Bottom panels, from left to right: cluster tidal radius R, global mass function slope a, central 
concentration c, and metallicity [Fe/H]. All these parameters are adopted from the on line catalog of Baumgardt & Hilker (2018), except log(tre), c and [Fe/H] are 
adopted from Harris (1996). The color-coded symbols, solid lines and texts have the same meanings as in Figure 3. For simplicity, we have omitted the confidence 


curves of the best-fitting functions. 


Figures 8(d)-(f), with vertical ranges indicate y-ray emissivity 
errors and the horizontal ranges roughly denote the extent of the 
boxy bulge (~10°) and nuclear bulge (~200 pc). We also 
displayed the best fit of GCE emission per unit stellar mass as 
red dashed horizontal line (є, = 1.8 x 10?’ erg s~'M5', Bar- 
tels et al. 2018) in the figures. 
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From Figures 8(d)-(f), it can be seen that GICs have a y-ray 
emissivity ~ 10-1000 times of higher than the Galactic boxy bulge 
and nuclear bulge, which indicates that tidally disrupted GICs may 
enhance the ^-ray emissivity of the Inner Galaxy greatly. More 
importantly, as GICs inspiral in toward the GC, it seems that 
normal stars are more likely to be stripped off the clusters by 
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Figure 6. Dependence of the y-ray emissivity on various physical properties of СІС. The layout of the panels, color-coded symbols, lines and comment texts in each 


panel are same as in Figure 5. 


Galactic tidal force, while MSPs are preferentially to be deposited 
to the Inner Galaxy, thus resulting in a negative dependence of c, 
On Reo, Rgc,per апа Кс аро in GICs. These features suggest that, in 
addition to the tidal stripping, mass seggregation also may plays an 
important role in allocating the objects that can be delivered to the 
GC. Therefore, an estimation of GIC contribution to the GCE must 
take these effects into account. We leave the implication of this 
problem to be addressed in Section 4.4. 

We end this section by emphasizing that the GIC parameters 
are not independent of each other. For example, there is a 
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significant correlation between Г and M in GICs (Cheng et al. 
2018), and massive GICs are usually found to have larger 
central velocity dispersion o,, escape velocity Vesce, core density 
Pc, and smaller half-mass (core) relaxation time 1, (tc). On the 
other hand, significant dependencce of cluster paramerters on 
К. and А. has been reported in literature (Djorgovski & 
Meylan 1994). These mutual relations suggest that the GIC 
parameters are highly coupled with each other, and they are 
closely connected with the dynamical evolution of GICs in the 
Milky Way tidal field (Meylan & Heggie 1997). In Table 3, we 


ҮІ 


ТаЫе 3 
Tested Correlations and Coefficients 


Correlations Total GCs Normal GCs Core-collapsed GCs 
GC param. L, E M Rgc L, €, M M 
L, 1.000 0.610 0.454 —0.267 1.000 0.544 0.464 0.259 
4.8e-5 0.004 0.105 es 0.004 0.017 0.417 
ey 0.610 1.000 —0.359 —0.559 0.544 1.000 —0.396 —0.343 
4.8e-5 © 0.027 2.бе-4 0.004 Ut 0.045 0.276 
M 0.454 —0.359 1.000 —0.193 0.464 —0.396 1.000 1.000 
0.004 0.027 es 0.015 0.017 0.045 © © 
Г 0.583 0.118 0.625 0.107 0.697 0.270 0.529 0.699 
1.5е-4 0.456 3.5e-5 0.529 1.1е-4 0.192 0.007 0.011 
A 0.454 0.368 0.141 —0.104 0.532 0.616 —0.067 0.517 
0.005 0.025 0.406 0.539 0.006 0.001 0.751 0.085 
Re —0.121 —0.255 —0.286 0.553 —0.390 0.621 —0.349 0.119 
0.475 0.128 2.5e-4 4e-14 0.054 9.3e-4 4.6e-5 0.538 
Dc 0.374 0.096 0.557 —0.525 0.658 0.352 0.651 0.152 
0.023 0.572 2e-14 le-12 3.5е-4 0.085 0 0.431 
B. —0.069 —0.535 —0.079 0.499 —0.275 —0.790 —0.202 0.581 
0.689 7.7е-4 0.349 Зе-10 0.194 4.3e-6 0.033 9.4e-4 
Oc 0.574 —0.217 0.911 —0.378 0.672 —0.189 0.942 0.763 
2.0e-4 0.197 0 9.0e-7 2.4е-4 0.367 0 1.5е-6 
Ry —0.363 —0.513 —0.226 0.666 —0.404 —0.548 —0.318 0.637 
0.027 0.001 0.004 0 0.045 0.005 2.2е-4 2.0е-4 
Ph 0.641 0.141 0.665 —0.559 0.728 0.124 0.767 —0.117 
2.0e-5 0.395 0 2e-14 3.6e-5 0.555 0 0.547 
teh —0.053 —0.609 0.271 0.559 —0.043 —0.657 0.245 0.797 
0.756 6.3e-5 5.5e-4 2e-14 0.839 3.6e-4 0.005 2.2e-7 
Vesc 0.559 —0.225 0.893 —0.389 0.652 —0.191 0.929 0.724 
3.3e-4 0.181 0 4.0e-7 4.1e-4 0.361 0 9.1e-6 
R, 0.039 —0.634 0.344 0.805 0.128 0.661 0.351 0.535 
0.818 2.5e-5 9.0e-6 0 0.543 3.2e-4 4.2e-5 0.003 
a 0.235 0.730 —0.286 —0.333 0.187 0.764 —0.279 —0.534 
0.161 2.9e-7 2.6e-4 1.8e-5 0.372 8.9e-6 0.001 0.001 
с 0.084 0.125 0.556 —0.103 0.267 0.177 0.607 0.206 
0.609 0.454 Зе-14 0.194 0.179 0.386 2е-14 0.283 
[Ее/Н] 0.482 0.611 —0.056 —0.439 0.494 0.586 —0.019 —0.163 
0.003 7Ле-5 0.497 2.1е-8 0.012 0.003 0.841 0.398 
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Table 3 
(Continued) 
Correlations Total GCs Normal GCs Core-collapsed GCs 
GC param. L, є, M Rgc L, €, M Rec L, €, M Кес 
X —0.425 —0.459 —0.068 0.680 —0.474 —0.490 —0.062 0.684 —0.392 —0.434 —0.140 0.839 
0.007 0.004 0.685 2.7е-6 0.012 0.011 0.764 1.2е-4 0.208 0.159 0.665 6.4e-4 
Y —0.169 0.129 —0.267 —0.198 —0.282 0.132 —0.430 —0.210 0.203 0.189 0.077 —0.238 
0.307 0.439 0.107 0.235 0.154 0.519 0.028 0.304 0.527 0.557 0.812 0.457 
Z —0.001 —0.107 0.121 —0.117 —0.126 —0.252 0.094 0.037 0.224 0.322 —0.301 —0.699 
0.998 0.523 0.470 0.483 0.530 0.214 0.648 0.858 0.484 0.308 0.342 0.011 
Rgc —0.267 —0.559 0.287 1.000 —0.302 —0.609 0.249 1.000 —0.497 —0.573 0.091 1.000 
0.105 2.6e-4 0.081 0.134 9.5e-4 0.220 0.101 0.051 0.779 
Крег —0.433 —0.562 0.143 0.662 —0.358 —0.419 0.052 0.663 —0.566 —0.944 0.504 0.685 
0.009 3.0e-4 0.400 7.9e-6 0.079 0.037 0.807 3.0e-4 0.055 3.9e-6 0.095 0.014 
Racapo —0.326 —0.604 0.280 0.974 —0.384 —0.666 0.237 0.966 —0.462 —0.622 0.203 0.951 
0.049 7.5e-5 0.093 0 0.058 2.8е-4 0.255 5e-15 0.131 0.031 0.527 2.0e-6 


Note. For each correlation, the Spearman's rank cofficient is shown in the upper panel of the grid, while the corresponding random correlation p-value is given in the bottom panel. 
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summarized the Spearman's rank coefficients of L,, €, M and 
К.с on various GIC parameters, the corresponding random 
correlation p-values are also listed in the table. 


4. Discussion 


By exploring the dependence of L, and €, on various GIC 
parameters in Section 3, we have established many correlations 
and trends between GIC dynamical evolution and their binary 
burning products. In particular, the meaured ^-ray emission is 
mainly contributed by the MSPs, which facilitates us to trace 
the dynamical formation of MSPs within GICs, their migration 
with host clusters in the Milky Way, and the final settlement 
and spatial distribution surrounding the GC. As introduced in 
Section 1, these processes are particularly important for the 
MSP interpretation of the GCE. Here we discuss our findings 
and conclusions as follows. 


4.1. GIC Dynamical Evolution and Formation of MSPs 


As self-gravitating systems, the dynamical evolution of GICs 
is balanced by the production of energy in the core and the 
outflow of energy from the cluster, and two-body relaxation is 
the fundamental process that dominates the transportation of 
energy and mass in GICs. Stars are driven to reach a state of 
energy equipartition by two-body relaxation, massive stars (or 
binaries) therefore tend to lose energy and drop to the cluster 
center (thus influences the types of binary burning products), 
whereas lower-mass stars tend to gain energy and move faster, 
and they will migrate outward (thus leading to the energy 
outflow in GICs) and drives the cluster envelope to expand 
(Heggie & Hut 2003). Binary burning is thought to be the 
internal energy source that supports the cluster against 
gravothermal collapse, and the energy production rate of GICs 
is sensitive to the encounter rate of the hard binaries (Cheng 
et al. 2018): 


2 
I, x f. A [20 “av, (1) 


where f, is the fraction of stars in binary, Ap e af c? is the 
binary encounter cross section, with a being the orbital 
separation of the binary. In practice, it is hard to determine 
the distribution of a by observation, while f, is found to vary 
not too much among GICs (i.e., fp ~ 196-2096; Milone et al. 
2012). Therefore, Г, is usually simplified as T œ Jp lo in 
literature. 

As binaries are disrupted or been driven to evolve to harder 
systems by dynamical ecounters, both f; and Aj will become 
smaller and the energy production rate of GICs may derease 
according to Equation (1). However, this may not happen in a 
real cluster. In fact, GICs contract their cores to increase the 
stellar density, which enhances the stellar encounter rates (i.e., 
Г and A) and thus the energy production rate. Accordingly, the 
cluster two-body relaxation time was adjusted to smaller 
values, to enhance the energy outflow rate and maintain the 
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balance between the energy production in the cores and the 
outflow of energy from the systems. As a result, the dynamical 
evolution of GICs will become more and more rapid, untill 
dynamical ejection of close binaries become important (see 
Section 4.2) and the cluster evolve into a core-collapsed GIC. 

With larger mass than normal stars, NS are expected to 
concentrate to cluster center and take part in the binary burning 
process, thus leading to the formation of LMXBs and MSPs in 
GICs. However, recent N-body simulations suggest that many 
GICs may host significant population of primordial stellar mass 
black holes (BHs), even evolve to an age larger than ~12 Gyr 
(Breen & Heggie 2013; Morscher et al. 2013, 2015; Wang et al. 
2016). Compared with NS, the BHs are expected to drop to 
cluster center first and form a high-density BH subsystem 
(BHS), where the BH burning process (ie., the dynamical 
hardening of BH binaries, Kremer et al. 2020b) may dominate 
the energy production of GICs and lead to the formation of 
gravitational wave sources (Rodriguez et al. 2016; Askar et al. 
2017; Anagnostou et al. 2020; Antonini & Gieles 2020). The 
thermal coupling of BHS with GICS is sufficient to reshape the 
structure of the cluster (Merritt et al. 2004; Mackey et al. 
2007, 2008; Giersz et al. 2019) and quench or slow down the 
mass segregation of less-massive objects (such as NS, binaries 
and BSS) in GICs (Fragione et al. 2018b; Weatherford et al. 
2018). Besides, through exchange encounters involving a BH, 
binaries are preferentially been transformed into BH binaries, 
BHS therefore may suppress the formation of other binary 
burning products in GICs, and clusters with a large number of 
BHs are more likely to have a lower formation efficiency of 
LMXBs, MSPs, CVs, ABs, and BSS (Fragione et al. 2018b; Ye 
et al. 2019; Kremer et al. 2020a). 

The mass segregation delay of heavy objects in GICs has 
been confirmed by the radial distribution study of X-ray 
sources, MSPs and BBS in 47 Tuc and Terzan 5 (Ferraro et al. 
2012; Cheng et al. 2019b, 2019a), where massive objects (i.e., 
MSPs, Bright X-ray sources or BSS) are found to drop to 
cluster center earlier and be more centrally concentrated than 
the less massive ones (ie., Faint X-ray sources and normal 
stars). The exception is M28, Cheng et al. (2020) reported an 
abnormal deficiency of X-ray sources in the cluster central 
region with respect to its outskirts, and argued that an early 
phase of primordial binary disruption by BHS may have 
happened in M28 (Cheng et al. 2020). On the other hand, 
significant evidence of BH burning has been reported in GIC w 
Cenaturi, in which the dynamical formation of X-ray sources is 
found to be highly suppressed, and BHS is the essential internal 
energy source that supports the cluster against collapse (Cheng 
et al. 2020). Indeed, with ~4.6% of the cluster mass being 
invisible, w Cenaturi is thought to be the cluster with the 
heaviest BHS in the Galaxy (Baumgardt et al. 2019). In 
Table 1, the y-ray emissivity of w Cenaturi is also the smallest 
among the 39 GICs, suggesting a low formation efficiency of 
MSPs in this cluster. 
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Figure 7. Dependence of GIC y-ray luminosity on Galacitic environment parameters. Top panels for cluster coordonates: with X point from GC toward the Sun (a), Y 
in direction of the Solar motion (b) and Z point toward the North Galactic pole (c). Bottom panels for the cluster orbital parameters in the Milky Way: 
Ry = yX ? 4 Y? + Z? represents the current GIC distance to the GC (d), while Rgc.per апа Rgc,apo Mark the perigalactic (e) and apogalactic (f) distances of the GIC 
orbit to the GC. The color-coded symbols and texts have the same meanings as in Figure 3. 


The evolution fate of the BHS is beyond the scope of this Lyx ee ag dynamically normal GICs (Figure 4(a)). 
paper. It is very likely that the majority of the BHs will be Besides, the GIC y-ray emissivity €}, a proxy of the MSP 


ejected from the cluster gradually, and host GICs tend to abundance of GICs, is highly correlated with the cluster 
contract their cores and increase the stellar density (Arca Sedda specific encounter rate A, with €, CX gc ean Figure 4(b). 
et al. 2018a). As a consequence, NS starts to concentrate to These correlations provide strong evidence for the dynamical 
cluster center and gradually take over the binary burning formation of MSPs in GICs. 

process, which then results in the formation a large number of On the other hand, by examing the dependence of L, and e, 
LMXBs and MSPs in GICs. Indeed, we have demonstrated that on the cluster structure paramters, we have tested the relation- 
L,, the proxy of the population of MSPs in GICs, is highly ship between MSP formation and the dynamcial evolution of 
correlated with the cluter stellar encounter rate, T, with GICs. For example, both L, and e, are anti-correlated with GIC 
7 According to Equation (1), the dynamical encounters between invisible yas еш Кова ap ш кышы к 
compact objects are nót included inT, Thus the energy input of BH burning are GICs (Figures 5 and 6). In particular, the ERe and 6у-Кһ 
also not counted here. correlations can be expressed as є x Ro 0934027 and 
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Figure 8. Dependence of GIC ^-ray esmissivity on Galacitic environment parameters. The layout of the panels are same as in Figure 7. The color-coded symbols and 
texts have the same meanings as in Figure 3. The measured ^-ray emissivity of the Galactic boxy bulge and nuclear bulge (Bartels et al. 2018) are shown as blue and 
olive rectangles in the botom panels, while the red dashed horizontal line represents the best fitting y-ray esmissivity of the Galactic boxy bulge and nuclear bulge. 
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in Figures 6(a) and (e), respectively. Further- 
more, Г is strongly dependent on the cluster stellar density р, 
with Го J p. /c. We have also found strong correlations 
between L, and р, and p, in dynamically normal GICs, with 
Ly о 003920. 14 апі L, x co in Figures 5(b) and (f). 
These trends are in agreement with the dynamical evolution of 
GICs, that is, as GICs contract their central regions to smaller 
radius and become more dense, more MSPs are expected to be 
dynamically formed in the clusters. 

Finally, it is necessary to emphasize that the dynamical 
formation of MSPs is far from finished in many dynamically 
normal GICs. As confirmed by the radial distribution studies of 
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X-ray sources (Cheng et al. 2019b, 2019a, 2020) and BSS 
(Ferraro et al. 2012) in GICs, it is possible that there are 
considerable numbers of primordial binaries and NS exist in the 
outskirts of the dynamically young GICs. Under the effect of 
mass segregation, these objects tend to drop to the dense core 
of the cluster, where they could be transformed into LMXBs 
and MSPs by the binary burning process, thereby increasing 
the abundance of MSPs in dynamically normal GICs. As a 
consequence, коше normal GICs with an older dyna- 
mical age (1.е., ty X tye lor ta X tin 5 are more likely to exhibit a 
larger abundance of MSPs, with the €, — һ correlation can be 
written as €, c 1,991591! in Figure 6(c), and the €,-th 
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correlation can be expressed as e, x 15082020 


respectively. 


in Figure 6(g), 


4.2. Dynamical Disruption of LMXBs and Ejection of 
MSPs in Core-Collapsed GICs 


From Equation (1), it can be seen that the evolution of GICs 
to core collapse is sensitively depends on the encounter cross 
section (i.e., Ap x a/ о?) of the binaries. With larger Аһ, the 
primordial binaries are expected to take part in the binary 
burning process first and are been transformed into close 
binaries, which may lead to the formation of many exotic 
objects, such as LMXBs (e.g., through exchange encounters 
involving NS), MSPs, CVs, ABs and BSS in GCs. However, 
when primordial binaries were exhuasted, GICs tend to contract 
their cores and increase the stellar density, thereby increasing 
the specific encounter rate (A) of stars in the cluster core. The 
enhancement of A is helpful for the extraction energy from 
harder binaries, and even very close binary systems, such as 
LMXBs, MSPs, CVs and ABs, may suffer strong dynamical 
encounters in core-collapsed GICs. Nevertheless, it is argured 
that the "burning" of very hard binaries is not sufficient to 
terminate the deep core collapse of GICs. The cluster cores are 
expected to contract further, untill dynamical interactions 
between single stars take effect and lead to the formation of 
new binaries in GICs (i.e., binary formation either via the “two- 
body binaries” or the “three-body binaries” channels, Heggie & 
Hut 2003). The sudden introduction of a large number of 
binaries in extremely dense core is sufficient to reverse the 
cluster core collapse to core bounce, and trigger gravothermal 
oscillations in GICs. 

Before the discussion of the dynamical evolution fate of 
LMXBs and MSPs in core-collapsed GICs, it is necessary to have 
an intuitive understanding on the dynamic effects of binary 
burning process of GICs. Accoding to the Hills-Heggie law, 
through binary-single encounters, hard binaries tend to increase 
their binding energy (E, = Gn? /2a) and become harder. The 
average increasement of E, per encounter is roughly proportional 
to its inital value, with ôE, 0.4-0.5 E, for equal-mass 
encounters (Heggie 1975; Hills 1975). The release of the binding 
energy is transformed into the kinetic energy of the ejected star, 
and the remaining binay will receive a recoil velocity as well, 
with vec x JOE, x 1 / Ja. Therefore, hard binaries tend to 
receive larger recoil velocity after strong encounters, which may 
lead to the recoil of the binaries out of the core into the halo, and 
sometimes out of the host cluster (Hut et al. 1992). 

With typical central escape velocity of 10km s! SZ 
Vese < 100 kms !, LMXBs and MSPs may receive a large 
recoil velocity and escape the core-collapsed GICs after 
suffering strong encounters. On the other hand, the specific 
encounter rate A of deep core collapse GICs is extremely high, 
it is possible for LMXBs to be dynamically disrupted (i.e., via 
the binary-binary encounters) or be greatly modified (i.e., 
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exchange encounters, etc.) before they otherwise could evolve 
into MSPs (Verbunt & Freire 2014). The net effect 1s that the 
population of MSPs tend to be gradually exhuasted in core- 
collapsed GICs. Indeed, compared with dynamically normal 
GICs, the core-collapsed GICs are found to exhibit smaller L, 
and є, with decreasing Re, trc, and increasing pe in Figures 5(a) 
—(c) and 6(a)-(c). Two well-known core-collapsed GCs, NGC 
6397 and NGC 6752, are proved to be the clusters with the 
smallest L, and e, in Table 1. 

The decreases of MSP populations in core-collapsed GICs 
have also been reported by de Menezes et al. (2019), which was 
interpreted as the disruption of LMXBs by the extremely large 
specific encounter rates in the dense cores. However, these 
authors have neglected the importance of MSP ejection in deep 
core collapse clusters. The lifetime of MSPs (~10'° yr) is at 
least one order of magnitudes larger than the cluster core 
relaxation time (tre ~ 105° yr), and once MSPs аге formed via 
the recycling process of LMXBs, it is hard to change their 
properties through direct stellar dynamical interactions. There- 
fore, the decline of L, and є, from dynamically normal GICs to 
core-collapsed GICs in Figures 5(a)-(c) and 6(a)-(c) may 
suggests an episodic ejection of MSPs during the deep core 
collapse phase, and the injection of MSPs into the Milky Way 
may plays an important role in shaping the spatial extent of the 
GCE. Nevertheless, it is also possible for some MSPs to be 
retained by the host cluster (Hut et al. 1992), and the retention 
of these MSPs may be responsible for the ~2 times larger 
detection rate of y-ray emission in core-collapsed GICs than in 
dynamically normal GICs by Fermi-LAT (Section 2.2). 


4.3. Tidal Stripping and Deposition of MSPs into the 
Galactic Center 


In addition to the internal effects, the structure of GICs is also 
subjected to the external influence of the host galaxy, which may 
affect the dynamical evolution of GICs in two ways. First, the 
expanding envelopes of GICs are expected to be truncated by the 
gravitational tidal field of the host galaxy, and compared with 
clusters in isolation, GICs in the tidal field may suffer an 
enhanced rate of energy outflow and mass loss. The tidal 
stripping therefore has a net effect in accelerating the dynamical 
evolution of the GICs. Second, through dynamical friction with 
background stars, GICs tend to lose energy and spiral in toward 
the center of the galaxy, where tidal stripping and interactions 
with the dense nucleus will eventually lead to the dissolution of 
the clusters (Tremaine et al. 1975). In particular, during the 
passages close to the Galactic bulge or passing through the dense 
disk, GICs may suffer gravitational shocks and these two effects 
could be enhanced greatly. The bulge shocking and disk 
shocking tend to heat up the GIC outskirts and increase 
evaporation of the cluster (Gnedin & Ostriker 1997), which in 
turn accelerate the core collapse and shorten the destruction time 
of GICs (Gnedin et al. 1999). Taking these aspects together, it is 
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suggested that GICs are the ancient building blocks of our 
Galaxy, and the inner galactic structures, such as the nuclear star 
cluster and the nuclear bulge, are assembled at least in part, by 
the merger of tidally disrupted GICs that spiraled into the deep 
gravitational well of the GC (Antonini et al. 2012; Antonini 2013; 
Gnedin et al. 2014). 

Observationally, the dynamical evolution of GICs and its 
coupling with host galaxy has been widely confirmed in 
literature, either via the stellar streams of tidal stripping (Leon 
et al. 2000; Sollima 2020), or the dependence of cluster 
structure parameters on GIC positions within the Galaxy 
(Djorgovski & Meylan 1994). As shown in Figure 9, with 
decreasing distance from the GC, GICs are found to be more 
concentrated and have smaller and denser cores. These trends 
are consistent with the dynamical evolution of GICs in the 
Galactic environment, namley, as GICs udergo orbital decay 
and spiral in toward the GC, the clusters tend to suffer stronger 
tidal stripping and will be more dynamically evolved. Besides, 
it is proposed that mass segregation and tidal stripping effects 
have a preference in evaporating the low-mass stars from GICs, 
which is strong enough to turn cluster initially increasing MF 
into MF that decrease toward the low-mass end (Vesperini & 
Heggie 1997; Baumgardt & Makino 2003; Lamers et al. 2013), 
and the global MF slope a is strongly anti-correlated with the 
half-mass relaxation time £j, in GICs (Baumgardt & Sol- 
lima 2017; Sollima & Baumgardt 2017). 

The relationship of GIC dynamical evolution and its 
coupling with the Galactic environment can also be confirmed 
by our y-ray data. For example, we showed that there is a mild 
anti-correlation between GIC ^-ray luminosity L, and their 
current distances Rg» perigalactic distances Rgc,per and 
apogalactic distances Rgcapo from the GC in Figure 7. 
Compared with the strong dependence of L, on Г in 
Figure 4(a), these correlations suggest that tidal stripping is a 
secondary factor affecting the binary burning process of GICs 
and the populations of MSPs formed therein. On the other 
hand, tidal stripping plays a crucial role in diminishing the 
mass of the GICs and leading to the final dissolution of the 
clusters, which can be further confirmed by the much stronger 
anti-correlations between GIC ^-ray emissivity є, and Ry, 
Recper ANd Кс аро in Figure 8. In fact, we also showed that there 
is no evident dependence of L, on GIC tidal radius R, in 
Figure 5(1), but є, is found to strongly anti-correlated with А, in 
Figure 6(1). Therefore, the loss of cluster mass is the most 
plausible reason for these anti-correlations. 

Finally, it is necessary to emphasize the importance of mass 
segregation in affecting the tidal stripping and dissolution of GICs. 
Unlike the low-mass normal stars that tend to migrate outward in 
GICs and will be stripped off the clusters preferentially, LMXBs 
and MSPs are more likely to segregate to GIC center and could be 
retained by the clusters for much longer time. As a result, GICs 
with increasing MS slope a are found to exhibit larger c, in 
Figure 6(j), and as GICs udergo orbital decay and spiral in toward 
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the Inner Galaxy, LMXBs and MSPs are more likely to be 
delivered to the GC, leading to a strong anti-correlations between 
€, and Rs, Rgc,per ANd Кс аро in Figure 8. 


4.4. Contribution of GIC MSPs to the GCE 


As introduced in Section 1, it is proposed that the GCE may 
arising from a large number (7107-107) of unresolved MSPs 
residing in the GC (Abazajian 2011; Yuan & Zhang 2014; 
Bartels et al. 2016). The generation of MSPs could be arise from 
the in situ star formation and evolution at GC (Yuan & 
Zhang 2014; Eckner et al. 2018; Gautam et al. 2022), or 
inherited from GICs that were brought to the Inner Galaxy by 
dynamical friction and assembled the Galactic nuclear star 
cluster and bulge (Bednarek & Sobczak 2013; Brandt & 
Kocsis 2015; Fragione et al. 2018a; Arca-Sedda et al. 2018b; 
Abbate et al. 2018). Compared with the dynamical channel 
within GICs, the in situ formation of MSPs at GC is disfavoured 
because of LMXBs, the progenitors of MSPs, are observed to be 
quite rare in the bulge of our Galaxy than the prediction of the 
MSP interpretation (Cholis et al. 2015; Haggard et al. 2017). 
Besides, Boodram & Heinke (2022) argued that the natal 
velocity kicks received by newly formed NS during supernova 
may lead to a lower abundance of MSPs (thus a lower synthetic 
y-ray surface brightness) in the central degrees (А < 150 pc) of 
the Galatic bulge with respect to its outskirts, and which is 
inconsistent with the measured ^-ray surface brightness profile of 
the GCE (Boodram & Heinke 2022). 

On the other hand, we have demonstrated that the dynamical 
channel of MSPs are very efficient during cluster dynamical 
evolution, and as GICs losing kinetic energy and sipiral in 
toward the GC, MSPs are preferentially to be deposited into the 
the deep gravitational potential well of the Galaxy,® which may 
avoid the problems reported in the work of Boodram & Heinke 
(2022). The spatial distribution of GICs are found to be 
spherically distributed around the GC (Baumgardt & 
Hilker 2018), and their velocity dispersion is highly isotropic 
within the central Galaxy (i.e., А 5 10 kpc; Vasiliev 2019). 
Therefore, the tidal stripping and ejection of MSPs from GICs 
to the Milky Way are also expected to create a spherically 
distributed MSP population in the Inner Galaxy, and the radial 
distribution profile resambles that of the GCE (Brandt & 
Kocsis 2015; Fragione et al. 2018a; Arca-Sedda et al. 2018b; 
Abbate et al. 2018). Moreover, considering the lifetime (~10'° 
yr) of MSPs is much longer than that (~10* yr) of LMXBs, the 
very small LMXBs to MSPs ratio observed in the Galatic bulge 
can be interpreted as the dynamical disruption of LMXBs 
during the deep core collapse phase of cluster evolution 
(Section 4.2), and the sudden interruption of the dynamical 
formation of LMXBs since host GICs were tidally disrupted 
(Brandt & Kocsis 2015). 


8 In fact, it can be seen from Figure 8(e) that the GIC MSPs could be 
effectively delivered to a distance of R — 100 pc from the Galactic center. 
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normal and core-collapsed GICs separately. The red, olive and blue texts donate the Spearman’s rank correlation coeficients of the core-collapsed, dynamically normal 


and the total GIC samples, respectively. 


To estimate of the contribution of Galactic GIC MSPs to the GICs, Hooper & Linden (2016) constrain the y-ray LF of MSPs 
GCE, the knowledge of the luminosity function (LF) of MSPs in GICs and found that they are as luminous as the Galactic 
is particularly important. By utilizing the cluster stellar field MSPs, with log-normal LF of Ly ~ 8.8 x 10? erg s! and 
encounter rate I' to estimate the formation rate of MSPs in о 0.62. Since MSPs would spin down quickly and fade away 
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in y-ray luminosity, Hooper & Linden (2016) therefore argued 
that the GIC MSPs can account for only a few percent or less of 
the observed GCE. Nevertheless, the LF of Hooper & Linden 
(2016) also predict that most of the GICs will be dominated by 
only one or two MSPs in y-ray luminosity, which is 
inconsistent with the radio survey result of MSPs in GICs, 
where many clusters detected by Fermi are found to host far 
more than one MSPs.” On the other hand, by using the cluster 
є; to scale the stellar mass deposited by GICs at GC, Brandt & 
Kocsis (2015) showed that a reasonable disrupted GIC mass, 
such as that calculated in Gnedin et al. (2014), may account for 
the total y-ray emission of the GCE. 

However, both the works of Hooper & Linden (2016) and 
Brandt & Kocsis (2015) have neglected the influence of cluster 
dynamical evolution on the formation rate of MSPs in GICs. As 
dicussed in Sections 4.1 and 4.2, the dynamical formation of 
MSPs is far from finished in dynamically young GICs, while the 
population of MSPs in core-collapsed GICs is underestimated, 
since strong encounters may have lead to the ejection of many 
MSPs from these systems. These effects can also be confirmed 
by the large variance of measured є, in GICs. From Figures 8(e)- 
(f), it can be seen that є, of GICs is 1—3 orders of magnitude 
larger than that of the Galactic bulge stars, the large variance of 
є; suggests that the final population of MSPs created by 
individual GICs could be more than tens times higher than that 
currently hosted in the cluster. Indeed, by assuming a cluster є, 
of log e, = 32.66 + 0.06 — (0.63 + 0.11)log M for GICs, 
Fragione et al. (20182) simulating the spirial in and tidal 
stripping of GICs in the Milky Way. Although their model has 
neglected the essential features of cluster dynamical evolution 
and stellar dynamics (i.e., internal energy sources, stellar mass 
segregation, core collapse, etc.) within GICs, the empirical 
є; — M relation of Fragione et al. (20182) is consistent with our 
results in Figure 3(b), and as GICs spiral in toward GC and get 
less massive, the c, of remaining cluster debris will become 
larger and larger, as displayed in Figures 8(e)-(f) More 
importantly, Fragione et al. (20182) showed that the cluster є, 
derived y-ray luminosity is about one order of magnitude larger 
than the observed GCE, and spin down of MSPs will reproduce 
a y-ray luminosity consistent with the GCE. 


5. Summary 


We present a y-ray study of the 157 Milky Way GICs based 
on the archival Fermi-LAT data with a time span of ~12 yr. By 
examing the dependence of GIC y-ray luminosity (L,) and 
emissivity (є, = L,/M) on various cluster parameters, our main 
findings are as follows. 

1. 39 GICs are found to be y-ray bright (i.e., with TS > 25) in 
our work, which corresponds to about 2596 (39/157) of the 
total Milky Way GICs. The detection rate of core-collapsed 


9 http:/ /www.naic.edu/ - pfreire/GCpsr.html 
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GICs (12/29) is ~2 times higher than the dynamically normal 
GICs (27/128). 

2. Compared with cluster mass (M), the GIC ^-ray emission 
is highly correlated with the stellar encounter rate (Г) and the 
specific encounter rate (A = Г/М), with pore and 
€, cx 978205 in dynamically normal GICs. These correlations 
provide strong evidence for the dynamical formation of MSPs 
in GICs. 

3. The formation of MSPs is also highly dependent on the 
dynamical evolution history of the host GICs. As GICs evolve 
to older dynamcial age (i.e., ty ос tg! or ty X tp") and become 
more compact (i.e., with smaller R, or А, and larger p, or pp), 
the dynamically normal GICs are expected to produce more 
MSPs, thereby exhibit larger L, and є, in Figures 5 and б. 

4. However, compared with dynamically normal GICs, the 
core-collapsed GICs are found to exhibit decreasing L, and є, 
as their cores undergo deep core collapse. This feature may 
imply that even LMXBs could be dynamically disrupted or be 
greatly modified in extremely dense cluster cores, and strong 
encounters may lead to the ejection of MSPs from core- 
collapsed GICs. 

5. With decreasing tidal radius (R,) and distance (К) from the 
GC, GICs are found to have increasing ^-ray emissivity, with 
Ey X g 0955 ang €, X дЕ Besides, є, is positively 
correlated with GIC stellar mass function slope o, with 
€, c 100957000 These correlations suggest that both tidal 
stripping and mass segregation effects are improtant factors 
influencing the abundance of MSPs in GICs, and as GICs 
undergo orbital decay and spiral in toward the GC, MSPs are 
more likely to be deposited into the GC rather than normal stars. 

6. We gauge the cluster €, is about 1—3 orders of magnitude 
larger than that of the Galactic bulge stars, which implies that 
GICs may enhance the ^-ray emissivity of the Galactic bulge 
greatly. The large variance of cluster e, may arise from the 
ongoing dynamical formation (or ejection) of MSPs in GICs, 
and the different degree of tidal stripping among the clusters. 
More importantly, the €, relations derived in our paper are in 
agreement with the empirical relation adopted in the simulation 
of Fragione et al. (2018a), which states that tidally disrupted 
GICs may provide a natural astrophysical explanation to the 
observed GCE. 


Acknowledgments 


This work is supported by the Youth Program of the 
National Natural Science Foundation of China No. 12 003 017. 


References 


Abazajian, K. N. 2011, Journal of Cosmology and Astro-Particle Physics, 
2011, 010 

Abazajian, K. N., & Kaplinghat, M. 2012, PhRvD, 86, 083511 

Abbate, F., Mastrobuono-Battisti, A., Colpi, M., et al. 2018, MNRAS, 473, 927 

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, Astronomy and 
Astrophysics, 524, A75 


Research in Astronomy and Astrophysics, 24:025001 (23pp), 2024 February 


Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33 

Abdollahi, S., Acero, F., Baldini, L., et al. 2022, ApJS, 260, 53 

Anagnostou, O., Trenti, M., & Melatos, A. 2020, PASA, 37, e044 

Antonini, F. 2013, ApJ, 763, 62 

Antonini, F., Capuzzo-Dolcetta, R., Mastrobuono-Battisti, A., et al. 2012, ApJ, 
31920, 111 

Antonini, F., & Gieles, M. 2020, MNRAS, 492, 2936 

Arca Sedda, M. 2020, CmPhy, 3, 43 

Arca Sedda, M., Askar, A., & Giersz, M. 2018a, MNRAS, 479, 4652 

Arca-Sedda, M., Kocsis, B., & Brandt, T. D. 2018b, MNRAS, 479, 900 

Askar, A., Szkudlarek, M., Gondek-Rositiska, D., et al. 2017, MNRAS, 464, L36 

Bahramian, A., Heinke, C. O., Sivakoff, G. R., et al. 2013, ApJ, 766, 136 

Bartels, R., Krishnamurthy, S., & Weniger, C. 2016, PhRvL, 116, 051102 

Bartels, R., Storm, E., Weniger, C., et al. 2018, NatAs, 2, 819 

Baumgardt, H., He, C., Sweet, S. M., et al. 2019, MNRAS, 488, 5340 

Baumgardt, H., & Hilker, M. 2018, MNRAS, 478, 1520 

Baumgardt, H., Hilker, M., Sollima, A., et al. 2019, MNRAS, 482, 5138 

Baumgardt, H., & Makino, J. 2003, MNRAS, 340, 227 

Baumgardt, H., & Sollima, S. 2017, MNRAS, 472, 744 

Bednarek, W., & Sobczak, T. 2013, MNRAS, 435, L14 

Belloni, D., Giersz, M., Askar, A., Leigh, N., & Hypki, A. 2016, MNRAS, 
462, 2950 

Belloni, D., Giersz, M., Rivera Sandoval, L. E., Askar, A., & CiecielA&Ag, P. 
2019, MNRAS, 483, 315 

Belloni, D., Zorotovic, M., Schreiber, M. R., et al. 2017, MNRAS, 468, 2429 

Boodram, O., & Heinke, C. O. 2022, MNRAS, 512, 4239 

Brandt, T. D., & Kocsis, B. 2015, ApJ, 812, 15 

Breen, P. G., & Heggie, D. C. 2013, MNRAS, 432, 2779 

Camilo, F., & Rasio, F. A. 2005, in ASP Conf. Ser. 328, Binary Radio Pulsars, 
328, ed. F. A. Rasio & I. H. Stairs (San Francisco, CA: ASP), 147 

Carlson, E., & Profumo, S. 2014, PhRvD, 90, 023015 

Chatterjee, S., Rasio, F. A., Sills, A., & Glebbeek, E. 2013, ApJ, 777, 106 

Cheng, K. S., Chernyshov, D. O., Dogiel, V. A., et al. 2010, ApJ, 723, 1219 

Cheng, Z., Li, Z., Fang, T., et al. 2019a, ApJ, 883, 90 

Cheng, Z., Li, Z., Li, X., Xu, X., & Fang, T. 2019b, ApJ, 876, 59 

Cheng, Z., Li, Z., Wang, W., et al. 2020, ApJ, 904, 198 

Cheng, Z., Li, Z., Xu, X., et al. 2018, ApJ, 858, 33 

Cheng, Z., Li, Z., Xu, X., et al. 2018, ApJ, 869, 52 

Cheng, Z., Mu, H., Li, Z., et al. 2020, ApJ, 892, 16 

Cholis, I., Evoli, C., Calore, F., et al. 2015, Journal of Cosmology and Astro- 
Particle Physics, 2015, 005 

Cholis, L, Hooper, D., & Linden, T. 2015, Journal of Cosmology and Astro- 
Particle Physics, 2015b, 043 

Clark, G. W. 1975, ApJL, 199, L143 

Clausen, D., Sigurdsson, S., & Chernoff, D. F. 2013, MNRAS, 428, 3618 

Daylan, T., Finkbeiner, D. P., Hooper, D., et al. 2016, PDU, 12, 1 

de Menezes, R., Cafardo, F., & Nemmen, R. 2019, MNRAS, 486, 851 

Di Mauro, M. 2021, PhRvD, 103, 063029 

Djorgovski, S., & Meylan, G. 1994, AJ, 108, 1292 

Eckner, C., Hou, X., Serpico, P. D., et al. 2018, ApJ, 862, 79 

Ferraro, F. R., Lanzoni, B., Dalessandro, E., et al. 2012, Natur, 492, 393 

Fragione, G., Antonini, F., & Gnedin, O. Y. 2018a, MNRAS, 475, 5313 

Fragione, G., Pavlík, V., & Banerjee, S. 2018b, MNRAS, 480, 4955 

Fregeau, J. M., Cheung, P., Portegies Zwart, S. F., et al. 2004, MNRAS, 352, I 

Fregeau, J. M., Gürkan, M. A., Joshi, K. J., et al. 2003, ApJ, 593, 772 

Freire, P. C. C., Abdo, A. A., Ajello, M., et al. 2011, Science, 334, 1107 

Gaggero, D., Taoso, M., Urbano, A., et al. 2015, Journal of Cosmology and 
Astro-Particle Physics, 2015, 056 

Gautam, A., Crocker, R. M., Ferrario, L., et al. 2022, NatAs, 6, 703 

Giersz, M., Askar, A., Wang, L., et al. 2019, MNRAS, 487, 2412 

Gnedin, O. Y., Lee, H. M., & Ostriker, J. P. 1999, ApJ, 522, 935 

Gnedin, O. Y., & Ostriker, J. P. 1997, ApJ, 474, 223 

Gnedin, O. Y., Ostriker, J. P., & Tremaine, S. 2014, ApJ, 785, 71 

Goodenough, L., & Hooper, D. 2009, arXiv:0910.2998 

Haggard, D., Heinke, C., Hooper, D., et al. 2017, Journal of Cosmology and 
Astro-Particle Physics, 2017, 056 

Harris, W. E. 1996, AJ, 112, 1487 

Heggie, D., & Hut, P. 2003, in The Gravitational Million-Body Problem: A 
Multidisciplinary Approach to Star Cluster Dynamics, ed. D. Heggie & 
P. Hut (Cambridge: Cambridge Univ. Press), 372 


23 


Feng et al. 


Heggie, D. C. 1975, MNRAS, 173, 729 

Heinke, C. O., Ivanov, M. G., Koch, E. W., et al. 2020, MNRAS, 492, 5684 

Hills, J. G. 1975, AJ, 80, 809 

Hooper, D., & Goodenough, L. 2011a, PhLB, 697, 412 

Hooper, D., & Linden, T. 2011b, PhRvD, 84, 123005 

Hooper, D., & Linden, T. 2016, Journal of Cosmology and Astro-Particle 
Physics, 2016, 018 

Hui, C. Y., Cheng, K. S., Wang, Y., et al. 2011, ApJ, 726, 100 

Hut, P., McMillan, S., & Romani, R. W. 1992, ApJ, 389, 527 

Ivanova, N., Heinke, C. O., Rasio, F. A., et al. 2006, MNRAS, 372, 1043 

Ivanova, N., Heinke, C. O., Rasio, Е. A., Belczynski, K., & Fregeau, J. M. 
2008, MNRAS, 386, 553 

Katz, J. I. 1975, Natur, 253, 698 

Kim, D.-W., Fabbiano, G., Ivanova, N., et al. 2013, ApJ, 764, 98 

Kremer, K., Chatterjee, S., Rodriguez, C. L., et al. 2018, ApJ, 852, 29 

Kremer, K., Rui, N. Z., Weatherford, N. C., et al. 2021, ApJ, 917, 28 

Kremer, K., Ye, C. S., Chatterjee, S., et al. 2020a, in IAU Symposium. 351, 
Star Clusters: From the Milky Way to the Early Universe, ed. A. Bragaglia 
et al. (Cambridge: Cambridge Univ. Press), 357 

Kremer, K., Ye, C. S., Rui, N. Z., et al. 2020b, ApJS, 247, 48 

Lamers, H. J. G. L. M., Baumgardt, H., & Gieles, M. 2013, MNRAS, 
433, 1378 

Leon, S., Meylan, G., & Combes, F. 2000, A&A, 359, 907 

Li, Z., Spitler, L. R., Jones, C., et al. 2010, ApJ, 721, 1368 

Lloyd, S. J., Chadwick, P. M., & Brown, A. M. 2018, MNRAS, 480, 4782 

Mackey, A. D., Wilkinson, M. L, Davies, M. B., et al. 2007, MNRAS, 
379, 140 

Mackey, A. D., Wilkinson, M. L, Davies, M. B., et al. 2008, MNRAS, 386, 65 

Merritt, D., Piatek, S., Portegies Zwart, S., et al. 2004, ApJL, 608, L25 

Meylan, G., & Heggie, D. C. 1997, A&ARv, 8, 1 

Milone, A. P., Piotto, G., Bedin, L. R., et al. 2012, A&A, 540, A16 

Moreno, E., Fernández-Trincado, J. G., Pérez-Villegas, 
Chaves-Velasquez, L., & Schuster, W. J. 2022, MNRAS, 510, 5945 

Morscher, M., Pattabiraman, B., Rodriguez, C., Rasio, F. A., & Umbreit, S. 
2015, ApJ, 800, 9 

Morscher, M., Umbreit, S., Farr, W. M., & Rasio, F. A. 2013, ApJL, 763, L15 

Petrović, J., Serpico, P. D., & Zaharijaš, G. 2014, Journal of Cosmology and 
Astro-Particle Physics, 2014, 052 

Pooley, D., & Hut, P. 2006, ApJL, 646, L143 

Pooley, D., Lewin, W. H. G., Anderson, S. F., et al. 2003, ApJL, 591, L131 

Ransom, S. M. 2008, 40 Years of Pulsars: Millisecond Pulsars, Magnetars and 
More, 983, 415 

Rasio, F. A., Pfahl, E. D., & Rappaport, S. 2000, ApJL, 532, L47 

Rodriguez, C. L., Chatterjee, S., & Rasio, F. A. 2016, PhRvD, 93, 084029 

Rodriguez, C. L., Morscher, M., Pattabiraman, B., et al. 2015, PhRvL, 115, 
051101 

Shara, M. M., & Hurley, J. R. 2006, ApJ, 646, 464 

Sivakoff, G. R., Jordán, A., Sarazin, C. L., et al. 2007, ApJ, 660, 1246 

Sollima, A. 2020, MNRAS, 495, 2222 

Sollima, A., & Baumgardt, H. 2017, MNRAS, 471, 3668 

Tam, P.-H. T., Hui, C. Y., & Kong, A. K. H. 2016, JASS, 33, 1 

Tam, P. H. T., Kong, A. K. H., Hui, C. Y., et al. 2011, ApJ, 729, 90 

Tremaine, S. D., Ostriker, J. P., & Spitzer, L. 1975, ApJ, 196, 407 

Vasiliev, E. 2019, MNRAS, 484, 2832 

Verbunt, F. 2003, in ASP Conf. Proc. 296, New Horizons in Globular Cluster 
Astronomy, ed. G. Piotto et al. (San Francisco, CA: ASP), 245 

Verbunt, F., & Freire, P. C. C. 2014, Astronomy and Astrophysics, 561, 
All 

Verbunt, F., & Hut, P. 1987, The Origin and Evolution of Neutron Stars, 
125, 187 

Vesperini, E., & Heggie, D. C. 1997, MNRAS, 289, 898 

Wang, L., Spurzem, R., Aarseth, S., et al. 2016, MNRAS, 458, 1450 

Weatherford, N. C., Chatterjee, S., Rodriguez, C. L., et al. 2018, ApJ, 
864, 13 

Wu, J. Н. K., Hui, C. Y., Wu, E. M. H., et al. 2013, ApJL, 765, 147 

Ye, C. S., Fong, W.-. fai., Kremer, K., et al. 2020, ApJL, 888, L10 

Ye, C. S., & Fragione, G. 2022, ApJ, 940, 162 

Ye, C. S., Kremer, K., Chatterjee, S., et al. 2019, ApJ, 877, 122 

Yuan, Q., & Zhang, B. 2014, JHEAp, 3, 1 

Zhang, P. F., Xin, Y. L., Fu, L., et al. 2016, MNRAS, 459, 99 


A., 


